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Abstract 

We present the program EvolFMC v. 2 that solves the evolution equations in QCD 
for the parton momentum distributions by means of the Monte Carlo technique 
based on the Markovian process. The program solves the DGLAP-type evolution 
as well as modified-DGLAP ones. In both cases the evolution can be performed in 
the LO or NLO approximation. The quarks are treated as massless. The overall 
technical precision of the code has been established at 5 x 10 -4 . This way, for the 
first time ever, we demonstrate that with the Monte Carlo method one can solve 
the evolution equations with precision comparable to the other numerical methods. 
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Programming language: C++ 
Computer: PC, Mac 
Operating system: Linux, Mac OS X 
RAM: less than 256 MB 
Number of processors used: 1 
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External routines /libraries: ROOT 
Subprograms used: 
Nature of problem: 

Solution of the QCD evolution equations for the parton momentum distributions of 
the DGLAP- and modified-DGLAP-type in the LO and NLO approximations. 
Solution method: 

Monte Carlo simulation of the Markovian process of a multiple emission of partons. 
Restrictions: 

(1) Limited to the case of massless partons. 

(2) Implemented in the LO and NLO approximations only. 

(3) Weighted events only. 
Unusual features: 

Modified-DGLAP evolutions included up to the NLO level. 

Additional comments: 

Technical precision established at 5 x 10 -4 . 

The EvolFMC version 1 was described in [1], but the actual code was not pub- 
lished. 

Running time: 

For the 10 6 events at 100 GeV: DGLAP NLO: 27s; C'-type modified DGLAP NLO: 
150s (MacBook Pro with Mac OS X v.10.5.5, 2.4 GHz Intel Core 2 Duo, gcc 4.2.4, 
single thread); 
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LONG WRITE-UP 



1 Introduction 

The evolution equations (EVEQs) for parton distribution functions (PDFs) 
and parton momentum distributions are one of the most efficient tools in 
calculating the radiative corrections in Quantum Chromodynamics (QCD) 
because they perform resummation of certain types of corrections up to infinite 
order. The PDFs are indispensable in any analysis of scattering processes 
which involve hadrons. The non-perturbative information on hadron structure, 
for the time being not calculable, is extracted from the experimental data 
and then used to create the PDFs at low energy scale. Among various types 
of evolution equations the most important is the family of the DGLAP-type 
equations [2]. Other widely used types of the EVEQs include BFKL [3], CCFM 
[4] or IREE [5]. In this paper we will discuss the DGLAP and modified-DGLAP 
types of EVEQs. However, a comment on the relation to CCFM EVEQs will 
be made. 

There are various numerical methods of solving the DGLAP-type EVEQs: 
Mellin transforms [6], evolution on a finite grid [7-9], expansion in Laguerre 
polynomials [10,11], expansion in Chebyshev polynomials [12,13], etc. The 
Monte Carlo (MC) methods differ from the other numerical techniques be- 
cause, in addition to providing the inclusive parton distributions, they supply 
also the complete tree of parton emissions during evolution. This allows one 
to construct the MC Parton Shower programs which provide the actual four- 
momenta of emitted quarks and gluons - a necessary input for any realistic 
analysis which must include experimental apparatus effects. One can find nu- 
merous implementations of the leading order DGLAP evolution in the MC 
Parton Shower codes; let us quote just a few examples: PYTHIA [14, 15], 
HERWIG [16,17], ARIADNE [18], GR@PPA [19,20]. So far the MC meth- 
ods have not been considered as a realistic alternative to the other numerical 
methods of solving EVEQs due to low precision and long time of computation. 
The presented here MC code EvolFMC is intended to fill-in this gap, profiting 
from the dramatic increase of the CPU power over two decades since NLO 
DGLAP evolution was formulated and solved numerically for the first time. It 
will be demonstrated in the following with the examples of numerical calcula- 
tions that EvolFMC can solve the (modified) DGLAP-type EVEQs with high 
precision (5 x 10~ 4 at least) within a reasonable CPU time. 

The program EvolFMC solves the EVEQs for the parton momentum distribu- 
tions by means of the MC simulation of the multiple emission of partons in 
the cascade. The emission process is of the Markovian type, i.e. each emission 
depends on the information from the previous emission only. The algorithms 
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are constructed on the basis of the Markovian process with simplified emission 
kernels which retain only the leading singularities. The complete kernels in the 
leading and next-to-leading approximation are then recovered by the standard 
reweighting procedure. 

The evolution is two-dimensional (x and t) by construction. However, the az- 
imuthal angle can always be added with a flat probability density distribution. 
Having identified the evolution time with certain kinematical variable, one can 
then reconstruct the four-momenta of all emitted partons. Such a procedure 
of reconstructing the four-momenta is, of course, exact only in the LO ap- 
proximation. In the NLO approximation the differential distributions will be 
strictly speaking correct only in the "inclusive" sense of the overall normalisa- 
tion. This is, however, the common (and in fact the only available) approach 
to the Parton Shower MCs. It should be mentioned here that recently there 
have been a few attempts to construct a true NLO Parton Shower algorithm 
with the help of "fully unintegrated" or "exclusive" partonic functions [21-24]. 

Apart from the standard DGLAP, two modified-DGLAP-type EVEQs can 
be solved by EvolFMC. The modifications involve a change of the argument 
of the coupling constant together with the introduction of a finite cut-off on 
its minimal value. The program works in the weighted mode only. The first 
version of the code, the EvolFMC v . 1, was described in [1], but the actual code 
was not published. The version v . 1 solved modified-DGLAP-type in the LO 
approximation only. The NLO evolution was available only in the standard 
DGLAP case. Also the structure of the code is rebuilt in version v. 2. 

Let us conclude this introduction with the following remark. The MC ap- 
proach to EVEQs has one important general advantage: with the help of the 
reweighting technique it is easy to introduce into evolution some additional 
effects or modifications. As an example let us mention the possibility of emu- 
lating the CCFM-type evolution. From the algorithmic point of view, having 
generated the azimuthal angles and reconstructing the transverse momenta 
of the partons in the shower, one can trivially construct the so-called "non- 
Sudakov" form factor and include it as a correcting weight. Of course, one 
has to remember that from the theoretical point of view the definition of the 
non-Sudakov form factor in the NLO case is a highly nontrivial problem; even 
the complete LO analysis is difficult in the context of the CCFM equation [25]. 
Additional modifications in evolution kernels can be done as well. 

The paper is organized as follows. In Section [2] we give a short theoretical 
overview of the evolution equations and their solutions by means of MC meth- 
ods. Section [3] describes in some detail the architecture of the presented MC 
code, EvolFMC version 2. Section [4] contains instructions on how to install 
the code on Linux and Mac OS X platforms. In Section [5] we present two 
demonstration programs included in the distributed version of the code. Sec- 
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tion [6] contains a detailed study of the technical precision of EvolFMC at the 
level of 5 x 10 -4 . A short summary in Section u\ concludes the paper. 



2 Theoretical background 



In this short theoretical overview we will only give a few formulae and defini- 
tions necessary to explain the notation, to define the problem and to present 
the solution. For detailed derivations and technical description of the algo- 
rithms we refer the reader to the extensive bibliography, which we will present 
here in more detail. The short and elementary description of the solution of the 
DGLAP-type EVEQs in terms of the Markovian process and its MC realiza- 
tion in the LO approximation has been presented in Ref. [26]. The extension to 
the NLO level and comprehensive description of MC algorithms that solve the 
Markovian evolutions of PDFs and parton momentum distributions has been 
given in Refs. [1,27]. The algorithms presented in [1] have been implemented 
in the first version of the code, EvolFMC v. 1. The extension to the modified- 
DGLAP-type evolutions and the detailed description of appropriate modified- 
DGLAP Markovian algorithms has been given in Refs. [28,29] for the LO case 
(implemented in EvolFMC v. 1) and in Ref. [30] for the NLO approximation. 
The implementation of algorithms from Ref. [30] as well as re-organization 
of the implementation of the DGLAP-type algorithms from Ref. [1] has been 
done in the second version of the code, EvolFMC v . 2, presented in this pa- 
per. Finally, in Ref. [13] the general and universal formalism of constructing 
Markovian algorithms, common for all DGLAP-type and modified-DGLAP- 
type evolutions, has been presented on the basis of the operator language. It 
is this formulation [13] that we will use in the rest of this section to describe 
the principles of the Markovian MC evolution. 

The evolution equation and its solution in the form of a master iterative 
formula for the Markovian MC algorithm can be expressed as follows: 



dtD(t) = K(t) D(t), i.e. d t D f (t,x) = ^ / dw X ff >{t,x,w)D f ,{t,w), 

/' 

(1) 
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and 



ED(t) 



dt\ / dt2 
t V Jti 



dU 



t-2 



... I dtx{ 'EK R (t N )G K v(tN, tN-i) + EG K v(t, t N _i)5t N =t \ 

JtN-l { J 



X 



x K R (t 2 )G K v (t 2 , h) + EG K v (t, ti)5 t2= 



x K fl (ti)G K v (ti, t ) + EG K y (t, t )^ 1= i ID (t ). 



Let us describe all the ingredients of eqs. Q and (§. 



(2) 



The multiplication of the matrices is understood as: J2f Jo ( ^' u; - 
Df(t,w) is the parton density function of the parton /. 
%ff(t,x,w) is the generalized evolution kernel built from the real and vir- 
tual parts: 



Xff/(t, x, w) =Kff/(t, x, w) + X R f,(t, x, w), 
Xj fl (t,x,w) = - 5 f f>5 x=w X v ff (t,x). 



(3) 



The operator E is defined as {E}f(x) = x, i.e. it turns parton distributions 

into parton momentum distributions, whereas summing and integrating over 

final degrees of freedom means in the MC language that we generate all 

possible final state configurations without any constraints. 

G K v is the solution of the evolution equation with the virtual kernel X^, (t, x. 

only 

{G K v{t,t')} ff {x,w) = 5 ff ,5 x=w e -WM. 



(4) 

$/(t, t'\x) is the Sudakov form factor, expressed in terms of the real emission 
part of the evolution kernel 



® f {t,t'\x) =J dt" X} f {t",x) = E/ dt " J — x ' Xf, f (t",x',x) 
t> f t> 

= J2^r f (t,t'\x). (5) 
/' 



The actual, normalized to unity, probability densities of the variables in each 
step of the Markovian process are visible in each of the lines of the eq. (|2]), 
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representing a single step in the emission chain: 

t 

H-l 

1 

_ g-^/i-iCMi-iNi-i) _)_ J d ^ e -®f i _ 1 (u,u-i\x i - 1 )^ 

e -* fi _ 1 <t.t i -i\'i-i) (6) 



In the program EvolFMC v . 2 we have implemented three types of the evolution 

'f'f 



differing by the definition of the evolution kernel 3C/y (t, x, w), X — A, B , C : 



X%P*\t,X,w) —6t+ln<f> x >ln\ 



aNLo(t + In <fix) 2z p R (°) 
2tt Z f ' f 



2zP^\z) 

(7) 



+ ( " jVLo(t 2 + ln ^ ) ) 2 2z(p / y(.) + Ap;f(.) 

where z = x/w. The parameter A is an arbitrary cut-off on the argument 
of the coupling constant, greater than Aqcd, necessary in order to avoid the 
singularity in the coupling constant. Note that the part of the real emission 
phase space excluded by the cut-off A is compensated for by the virtual form 
factor defined in eq. ^ as the integral over phase space of a real emission. As 
a consequence the momentum sum rule is preserved. ln<f)x takes one of the 
following three forms: 
A:ln0 x = O, (the DGLAP evolution) 

B': ln<p x = ln(l - z), (the modified-DGLAP B'-type evolution) 

C: ln0 x = ln(iu(l - z)), (the modified-DGLAP C'-type evolution). 

The term APpf X (z) is added to remove the double counting caused by the 
change of the argument of the coupling constant, according to the prescription 
of Ref. [31]. Namely, from the expansion 

OiNLo{t +ln</>) = 0>NLO 

one obtains 



APpy )D (z) = AP^ (z) = A>ln(l - z)p; 
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z 



Note that in the case C we use the counter term identical as in the B' case. 
The additional piece related to Inw is of a genuine beyond-DGLAP origin, 
i.e. it is absent in the DGLAP kernel. Therefore, there is no double counting 
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and no need to subtract it. The universal LO part P^f\z) and the NLO part 

Pf,\ (z) are given in [32,33]. Finally, the coupling constant at the NLO level 
has the standard form 

m 2tt / /3i ln(2t - 2 In A ) \ 

"loW = -5-77 — i — rr, a NLO {t) = a LO {t) 1 - «loW 



/3 (t-lnA )' JV1 ^ W 47r A 



(9) 



On the technical side, in the actual MC algorithms implemented in EvolFMC 
v. 2 we do not use the complicated kernels ([7]). Instead, a series of simplified 
kernels is introduced. Each of them is chosen in such a way that it 

retains only the leading singularities of the exact kernel % R<yX ^ while all the 
complicated but finite structure is temporarily discarded: 

xXpf\t,x,w)= a ^2zP?f\z), (10) 



for the X = A case and 

XAj,j (t,X,W) = — ZZ ^f'f l 2 J fc, t+ln</. x >lnA, (H) 

f JVL( 

77 



-Wl-V)^ „ 0\7.()(/ • lllOy)., - • 

2tt 



for the cases X = B' and X = C". The variable z = x/w. The constant M x ^° 
is defined as 



0, if P5 f 0) (^) ^0, 
/J U ifP / T(,) = 0, 

and 77 is a dummy technical parameter. The functions Af}(z) and FpJ(z) are 
a convenient parametrization of the full LO kernels Ppp{z) 

zP?f\z)=^8 flf Af} + Fp}(z), (13) 

see Appendix C of Ref. [1] for the complete list of them (for example, for the 
zP$°\z) kernel we have A$ = 2C F and F$(z) = C F {-2 - z - z 2 )). The 
NLO kernels zP R ^ (z) in the form used in the code are explicitly given in 
Appendix A of Ref. [1]. The exact kernel % R ( X ' is recovered at the end by the 
standard reweighting procedure. The correcting weight is 

,(n) =p $f n (t,t n \x n )-§f n (t,t n \x n ) 



X 



1=1 ^-/i/i-l x «' / 
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The product runs over all generated partons in a given MC event with multi- 
plicity n, and the form factor is constructed from 3Cff Xi, a^-i), 
in analogy to eq. ([5|. 

The actual expressions for the form factors $ and $ are fairly complicated, 
especially for the cases B' and C, and we will not quote them here, referring 
the interested reader to the original papers. Let us only remark that, as seen 
in eq. (|5]), the form factors are defined as two-dimensional integrals. In the 
case of the simplified form factor $ both integrals can be done analyticaltyp"] 
On the contrary, in the full form factor $ only one integration can be done 
analytically. The other one has to be done numerically on an event-per-event 
basis. 



3 Overview of the software structure 



The program EvolFMC is written in the C++ language. To compile and link the 
code we use the autotools utility. It allows us to compile/link the code on 
many platforms in a simple way. The code has been routinely compiled and run 
under the Linux and Mac OS X 10.5 operating systems. From the user point 
of view, the only difference between these two systems is in the compiler's 
options inside the file configure . in, which is included in the main folder 
of the project. The central part of the EvolFMC source code is the MarkovMC 
library located in the MarkovMC folder. This folder includes the essential source 
code necessary to solve evolution equations. This part of the code requires only 
basic C++ libraries and an external random number generator (RNG). We use 
the generator TRandom3 from the ROOT package as a default RNG. With this 
generator we have reached the precision below 0.05%. In case when ROOT is not 
available on a given system platform, one should replace in the wrapper class 
rndm the name TRandom3 with the name of the other RNG. A simple main 
program in the DemoO folder uses only the standard C++ libraries and can be 
built and executed without ROOT. On the other hand, the more sophisticated 
source code in the Demol folder of the distribution version of the project uses 
the ROOT library, mainly for booking, filling and drawing histograms and more, 
see below. 

The library in the MarkovMC folder has a modular structure in the sense that 
the algorithms that solve a particular type of the evolution equations are 
implemented in separate classes and located in separate source files. Each 



1 In fact this analytical integrability is one of the criteria in choosing the form of the 
simplified kernels (11 ). This is for example why the 1/(1 — z) term multiplies artifi- 
cially also the F- function in eq. (11). The constant M is added to avoid potentially 
dangerous zeroes at the NLO level. 
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class implementing the Markovian MC algorithm for a given type of evolution 
equation includes all formulae, in particular the Sudakov form factors, specific 
to a given evolution type. All classes specific to one type of evolution inherit 
from the common base class MarkovianGen. This modular structure also makes 
it easier to add in the future any new type of the QCD evolution of the parton 
distributions to the code. 

3.1 Structure of folders 

In the following the structure of the folders of EvolFMC is described: 

• MarkovMC - the folder containing the library of the Markovian MC engines. 
Source codes of the classes solving various types of the evolution equations 
are placed in separate files. 

• DemoO - the folder with the simple demonstration program DemoO written 
in the C language. This program demonstrates the standalone usage of the 
library MarkovMC, without the use of ROOT. 

• Demol - the folder hosting the demonstration program Demol, a template 
program for the advanced user of EvolFMC. Demol requires ROOT to be in- 
stalled in the system. The subfolder Demol /work contains scripts necessary 
to run Demol. 

• m4 - the folder containing a script which defines properly a path for the 
ROOT libraries. This script is used by the automake program. 

Each of the above folders contains also scripts (in the files Makefile . am), 
which are required by the automake utility. 

3.2 Source code in folder MarkovMC 

The source code of the MarkovMC library consists of three categories of files: 
(1) the files containing base classes of the Markovian MC generator, (2) the 
files which contain classes specific to a particular type of the QCD evolution 
equations, and (3) the files with some auxiliary classes. 

(1) Base classes 

• markoviangen. cxx, markoviangen.h: 

The class MarkovianGen contains the essential part of the Markovian 
MC algorithm, member functions and data members common to all 
types of the QCD evolution. In particular, it executes the main Marko- 
vian loop over parton emissions. Virtual member functions encapsulate 
evolution details. 

• kernels, cxx, kernels, h: This class defines the DGLAP LO kernels. 
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• kernels_nlo . cxx , kernels_nlo . h: This class defines the DGLAP NLO 
kernels; it inherits from the simpler class kernels. 

(2) Auxiliary classes 

• gaussintegral . cxx, gaussintegral .h: the standard Gauss integra- 
tion procedure, translated from the Fortran GNU library to C++. 

• rndm . cxx , rndm . h: the wrapper class of the random number generator; 
it is derived from the ROOT class TRandom3. 

(3) Classes implementing one particular type of the QCD evolution equation 
for the parton momentum distributions: 

• dglap_lo . cxx , dglap_lo . h: implementation of the DGLAP LO evolu- 
tion; this class is derived from the classes MarkovianGen and kernels. 

• dglap_nlo . cxx dglap_nlo . h: implementation of the DGLAP NLO evo- 
lution; this class is derived from the classes MarkovianGen and kernels_nlo. 

• bprim_lo . cxx, bprim_lo.h: implementation of the modified-DGLAP 
evolution scheme B', the LO case; this class is derived from the classes 
MarkovianGen, kernels and gaussintegral. 

• bprim_nlo . cxx , bprim_nlo . h: implementation of the modified-DGLAP 
evolution scheme B', the NLO case, the basic algorithm; this class is de- 
rived from the classes MarkovianGen, kernels_nlo and gaussintegral. 

• bprim_nlo_aux . cxx , bprim_nlo_aux . h: implementation of the modified- 
DGLAP evolution scheme B', the NLO case, the auxiliary algorithm 
(for tests only!); this class is derived from the classes MarkovianGen, 
kernels_nlo and gauslntegral. 

• cprim_lo . cxx, cprim_lo.h: implementation of the modified-DGLAP 
evolution scheme C, the LO case; this class is derived from the classes 
MarkovianGen, kernels and gaussintegral. 

• cprim_nlo . cxx , cprim_nlo . h: implementation of the modified-DGLAP 
evolution scheme C, the NLO case, the basic algorithm; this class is de- 
rived from the classes MarkovianGen, kernels _nlo and Gaussintegral. 

• cprim_nlo_aux . cxx , cprim_nlo_aux . h: implementation of the modified- 
DGLAP evolution scheme C, the NLO case, the auxiliary algorithm 
(for tests only!); this class is derived from the classes MarkovianGen, 
kernels_nlo and gaussintegral. 

3. 3 Inheritance pattern of classes of MarkovMC library. 

As already indicated, all the classes that are used to solve the evolution equa- 
tions for the parton momentum distributions by means of the Markovian MC 
method are derived from a few base classes. Depending on the type of the QCD 
evolution, the base classes are: MarkovianGen, kernels and its derived class 
kernels_nlo, and the auxiliary class gaussintegral. The derived classes im- 
plement details of the particular type of the QCD evolution equations (e.g. 
dglap_lo . cxx). In particular, these classes include member functions which 
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generate randomly the evolution time, the parton flavor and the z- variable, 
which are defined as virtual member functions in the base class MarkovianGen. 
The base class MarkovianGen implements all essential parts of the Markovian 
algorithm and a few auxiliary functions. The central member function of this 
class is GenerateEvent. The class MarkovianGen does not know the details 
of the evolution kernels - they are implemented in the class kernels and/or 
kernels_nlo. 

The gausslntegral class owns integration methods. These methods are nec- 
essary to calculate the Sudakov form factors for more complicated types of 
the QCD evolution. 

Depending on the complexity of the equation, the structure of inheritance has 
different forms. As an example we present in Fig. [I] the inheritance scheme for 
the most complicated case of the modified-DGLAP NLO C'-type evolution. 



Class MarkovianGen 



GenerateEvent 

GenerateZ 

GenerateT 

KernelWeight 



«EventGenerator» 



Class gausslntegral 



Integrate 



~ ~ ^ «lntegral» 



Class kernels 



GenerateFlavor 
LO Kernels 



Class kernels nlo 



NLO kernels 
Probability matrices 
T 



«NLO kernels» 



Class cprim nlo 



Fig. 1. The structure of the derivations in the case of the modified-DGLAP NLO 
C'-type evolution. 



3.4 General design o/MarkovMC library 



The classes of the MarkovMC library and the organization of the source code 
have been designed in such a way that (a) an infrastructure is available for any 
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kind of the Markovian MC implementing any kind of the QCD evolution, (b) 
it is easy to include or exclude any group of classes implementing any type of 
the QCD evolution. It is therefore not surprising that the classes which solve 
different evolution types are completely independent from each other. One can 
exclude a particular class from the code without loosing functionality of other 
classes. Also, adding more evolution types would not modify the structure of 
the library. In the practical application, if one needs, for instance, a solution of 
the modified-DGLAP NLO scheme C, then one should include only the header 
file cprim_nlo .h. A simple example will be shown in the DemoO program in 
the following sections. 



3.5 Description of base class MarkovianGen 



Virtual member functions. 

The member functions from this group are implemented in the derived classes: 

• void MarkovianGen: : GenerateEvent (double Set, double &Vx, double feweight, 
double tmax, int &Flavor, double epsTSolver = . 0001) - generates 

a single MC event according to the Markovian algorithm. 

• void MarkovianGen: : GenerateEvent () - generates a single MC event ac- 



cording to the Markovian algorithm, a "wrapper" function, see Sect. 5.2 for 
explanation. 

• double GenerateT (double rndm, double t_prev, double TStop, double 
epsTSolver, double Vx) - generates the evolution time^j 

• double GenerateZ (double rndm, double t, double TO, double epsTSolver, 
double Vx) - generates the actual light-cone variable z. 

• int GenerateFlavor (double rndm, int oldFlavor) - generates the parton- 
flavor index. 

• double KernelWeight (double t, double z, double Vx) - provides part 
of the MC weight turning the simplified kernel into the exact kernel. 

• double DeltaRealPart (double t_new, double t_old, double z)-part 
of the MC weight from the Sudakov form factor evaluated analytically. 

• double DeltaVirtualPart (double t_new, double t_old, double z) - 
part of the MC weight due to the Sudakov form factor evaluated numerically. 

• virtual void init() - initialization of an object. 

• void AddParticle (double t, double x, double z, int f , double weight, 
int index) - stores data of a single parton emission, a private function, 

not to be used outside the MarkovianGen class. 

• Other auxiliary member functions used to transmit the flavor type between 
the class MarkovianGen and the class kernels: 

2 The technical parameter epsTSolver sets a precision for the TSolver function 
used for some evolution types. 
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• void SetActualFlavor (int flavor) - transmits the actual parton fla- 
vor to kernels. 

• void SetOldFlavor (int flavor) - transmits the old flavor to kernels. 

• "Getter" functions for the external users: 

■ bool GetParticle (double &t, double &x, double &z, int &f, double 
feweight, int index) - provides the external user with the information 
about partons generated in the last MC event; index runs from to 
EventMultiplicity, the variables t, x, z, f and weight describe the emis- 
sion of the number index, weight is the cumulative weight; index = 
returns information on the initial parameters of GenerateEvent and in 
this case z = 1. 

■ int GetEventMultiplicityO - returns a number of particles generated 
in the last MC event. 

All parameters of the virtual functions belong to the following list: 

• double rndm - the random number, 

• double t - the evolution time, 

• double TO - the start of the evolution time, 

• double TStop - the end of the evolution time, 

• double z - z — x ncw /x i d , 

• double Vx - as input: the light-cone x id variable before the emission; as 
output: the x new variable after the emission, 

• double t_new - t new the generated current evolution time, 

• double t_old - t \d the evolution time of the previous emission, 

• double knew - / new the generated current flavor, 

• double kold - / id the flavor of the previous emission, 

• double flavor - the flavor (depends on the function - explanation above), 

• double weight - as input: the initial weight to be assigned to the event 
(for example from the generation of the initial condition); as output: the 
cumulative weight of the event, 

• epsTSolver - the precision for the TSolver function. 

3.6 Small parameters in classes o/MarkovMC library 

Constructors and other methods of the classes in the MarkovMC library have as 
formal parameters several small parameters, which we call epsilon-parameters. 
They may be of technical or physical character. Let us explain them with 
explicit examples. 

(1) The constructor 

dglap_nlo (double TO, double lambda, double NumberOf Flavor , 
rndm *rnGen, double epsilon_IRC = 0.0001, 
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double epsilon_Zmin = 0.0001) 



contains two physical epsilon-parameters: 



epsilon_IRC is the dummy infrared cut-off at z = 1 for the DGLAP evo- 
lution. Solutions of the evolution equations do not depend on its value as 
long as it is kept small enough (10~ 4 by default). 

epsilon_Zmin is the minimal value of the final x- variable. In practice it is 
used as a cut-off for the z-variable in the kernels which exhibit a logarithmic 
divergency in the small z limit. Note that the MC program will generate 



the distribution for x < e Zmin but it will be incorrect, see Sect. |5 . 1 1 for more 
comments. 



(2) In more advanced classes, such as cprim_nlo, there is another technical 
epsilon-parameter: epsTSolver. It is used in the method GenerateEvent of 
the base class to set the precision of the important TSolver member function 
which inverts numerically an arbitrary one-dimensional function. 

(3) The last two technical epsilon-parameters are defined in the auxiliary class 
gausslntegral, see for instance one of its member functions 

double dqags (double a, double b, double epsabs, double epsrel, 
double feabserr, int feneval, int &ier) , 

where epsabs and epsrel are used to set the technical absolute and relative 
precision of the numerical integration. 



3. 7 Initial parton momentum distributions 



The MarkovMC library does not generate the initial parton momentum distri- 
butions - the user is supposed to provide the initial values of x- and flavor- 
variables for the GenerateEvent method. The actual generation of the initial 
parton momentum distributions is therefore done by an external MC appli- 
cation. We provide two examples of such an external environment. The first 



one, in the folder DemoO (see Section 5.1), simply uses fixed values of starting 



Xstart and Flavor. The second, a more advanced example in the folder Demol 



(see Section 5.2), uses the adaptive MC generator TFoam (part of the ROOT sys- 
tem) to generate the initial densities. In Demol we use the gluon (G) and quark 
singlet (Q) PDFs with three massless quarks in the following notation: 



Dq 



E(a, + a,) (is) 



15 



and 






Dl{x) = Dl(x) = D° s (x) = ^D° sea (x), 
D° sea (x) + DlJx) + Dl ai (x). 



(16) 



More details on the actual implementation of the TFoam-based generation can 
be found at the end of Section 15.21 



4 Installation instructions 

This section instructs the user of EvolFMC v. 2 how to install the program on 
two system platforms: Linux and Mac OS X. The syntax of Linux commands 
is given for the bash shell. The first few steps concern installation of the ROOT 
package. As explained earlier, the library MarkovMC as such does not need 
ROOT. However, in the distributed version the ROOT package is required as a 
source of the random number generator for the library MarkovMC. 

The more advanced demonstration program Demol exploits ROOT as a his- 
togramming package and uses its persistency mechnism. 

A step-by-step installation procedure of EvolFMC looks as follows: 

(1) Check if ROOT is installed in the system and find its location. In the case 
there are several versions of ROOT in the system, choose the preferred one. 

(2) Check if the environmental variable ROOTSYS is defined correctly: echo 
$R00TSYS. If several versions of ROOT are in the system, define the ROOTSYS 
variable path to the correct /preferred version: 

export ROOTSYS=path_to_your_root 

Check if the shell variable LD_LIBRARY_PATH (or DYLD _L I BRARY _P ATH in 
Mac OS X) contains a correct path to the ROOT's library. If not, then 
execute: 
under Linux: 

export LD _L I BRARY _P ATH= $ LD _L I BRARY_P ATH : $R00TSYS/lib 
under Mac OS X: 

export DYLD_LIBRARY_PATH=$DYLD_LIBRARY_PATH : $R00TSYS/lib 

It is convenient to put this command into the bash-shell configuration 
files: .bashrc or . bash_prof ile. Finally, check if the $PATH variable 
includes a correct path to the ROOT binaries. 

(3) Add the path to the project to the variable LD .LIBRARY _P ATH (DYLD_LIBRARY_PATH 
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in Mac OS X). This is done under Linux by: 

export LD_LIBRARY_PATH=$LD .LIBRARY _PATH : path_to_the_pro j ect/lib 
or under Mac OS X : 

export DYLD_LIBRARY_PATH=$DYLD_LIBRARY_PATH : path_to_the_pro j ect/lib 

Note: For the demonstration programs DemoO and Demol the user may 
skip the above commands because this path is already set in the appro- 
priate Makefiles. 

(4) Build the program from the commad line in the main project folder: 
autoreconf -i — force 

under Linux: ./configure 

under MacOSX: ./configure — enable-platf orm=macos 
make 

(5) Test the correctness of the installation: 
(cd DemoO; . /verif y_benchmarks) . 

For more details on the above test as well as on how to run two demon- 
stration programs DemoO/Demo and Demol/DemolPr see the next section. 

The authors have also managed EvolFMC using two popular integrated software 
development packages: Kdevelop and Eclipse. Let us hint on how to initialize 
EvolFMC as a project within these development tools: 

(1) Kdevelop 

From the main folder of the project just type in the shell: 
kdevelopfe 

then from the menu <pro j ect> choose <import> and set <pro j ect type> 
to <Generic C++ Application (Automake-based)>. Next time you open 
Kdevelop, the configuration files will be already in place. 

(Do not forget in the menu <project>^<project options>^<conf igure options>: 
in the window <Conf iguration> to choose default instead of t/efcutpH) 

(2) Eclipse 

The configuration files are included, so it is enough to invoke Eclipse and 
from the menu <File>^<Import>^<General> choose <Existing pro j ects>. 
The list of existing projects should appear, including the current EvolFMC. 

4-1 Testing correctness of installation 

Finally, let us explain how to test quickly the correctness of the installation. 
This is done by means of executing a special "benchmark test" in form of 
the bash script verif y_benchmarks included in the subfolder DemoO of the 
distribution folder. This test compiles and links the program DemoO/demo . cxx, 
and then runs it in a sequence for all eight implemented types of the evolution. 
Text outputs from these runs, containing a printout of the variable x, the flavor 

3 Unless you really need debugging. 
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type and the MC weight for 100 events, are produced and compared using the 
dif f utility against the benchmark outputs stored by the authors of the code 
in the subfolder DemoO/bmarks .outputs. If the installation procedure and all 
the settings are correct, then there should be no differences between the stored 
and the current output disk files. The stored outputs have been generated on 
the system iMac Core 2 duo with gcc 4 . 2 under Ubuntu 8. 
Summarizing: 

(1) The benchmark test can be invoked as follows: 
cd DemoO 

with the help of the bash script 
. / verif y -benchmarks 

(2) If needed, the user may create his/her own new set of benchmark output 
files with the help of the bash-shell script in the DemoO folder: 

. /create_benchmarks 

Alternatively, the above benchmark can be executed with the command make 
bmark. 



5 Two demonstration programs 



In the following we describe in a more detail two demonstration programs in 
the subfolders DemoO and Demol, which the user should run after installation 
of EvolFMC. They are also meant as the templates for applications which use 
EvolFMC in studies related to the perturbative QCD - most likely as a testing 
tool for other MC programs implementing the QCD evolution of the parton 
momentum distributions, or as part of some bigger MC application. These two 
programs are already built during installation, see the previous section, and 
can be executed as follows: 

(1) A simple demonstration program: 
cd DemoO 

make start 

(2) A more advanced demonstration program: 
cd Demol/work 

make start 

Four histograms are recorded in the disk file. To visualize them: 
make plot 

Let us describe these two demo programs in a more detail. 
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5. 1 Simple demonstration program demo . cxx 



The simple demo . cxx program is included in the subfolder DemoO in order to 
demonstrate how to generate MC events for any of the eight evolution types 
supported by EvolFMC. The listing of the whole program demo. cxx is given 
in the Appendix. Let us present and explain the crucial instructions in the 
demo . cxx source code, using the example of just one evolution type - LO 
DGLAP: 

[. • J 

#include "rndm.h" (1) 

#include "dglap_lo.h" 

[...] 

int mainO 
{ 

[. . J 

rndm * RNgen = new rndm() ; (2) 
dglap_lo * dglap_ll = new dglap_lo (TO, Lambda, 

numberOf Flavors, RNgen, eps IRC, epsZmin) ; (3) 
for(int i=0; KnumberOf Events ; i++) 
{ 

dglap_ll->GenerateEvent (Tc , Xstart , weight , TStop , Flavor) ; (4) 
[. . J 

} 

[. • J 

> 

(1) The appropriate header files are included: one for the class of the random 
number generator and another one for the class of the chosen evolution 
type. 

(2) The object of the RNgen class being the random number generator is 
created. 

(3) Next, the MC generator object dglap_ll of the class dglap_lo is created. 
It is the central object of the above code and it is used to generate 
a series of the MC events. The arguments of the dglap_lo constructor 
are: the initial (starting) value of the evolution time Tmin, the value 
of A.qcd Lambda, the number of active flavors numberOf Flavors and the 
pointer to the object of the random number generator RNgen. In addition, 
the parameter eps IRC is the technical cut-off used for regularizing the 
distribution 1/(1 — x) + at x = 1 in the kernel JCffi of eq. (j3j). The final 
result will be independent of epsIRC, if it is kept small enough. The last 
parameter, epsZmin, is the minimal requested value of the x variable 
to be generated. A non-zero value of the cut-off e min is necessary only 
in the NLO cases, due to the presence of the (\nz)/z singularity in the 
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NLO DGLAP kernels. 1/z in the kernels is cancelled for evolution of the 
momentum distributions but In z is still present, hence some form of a cut- 
off is needed. Note that the solution for x > e min given by the program 



is always independent of this cut-off, see Sect. 3J3 for more details^] 
Note that both epsIRC and epsZmin have default values assigned by 
the constructor and their redefinition, as done in the above example, is 
optional. 

(4) Finally, the MC events are generated with the method GenerateEvent. 
The meaning of the arguments of GenerateEvent is the following: 
Tc is a parameter for technical tests, not to be used. The initial conditions 
of the evolution are set by Xstart and Flavor, being the initial values of 
the x- variable and the flavor typd^~| TStop is the maximal value of the 
evolution time, weight is the initial weight assigned to the event, nor- 
mally set to 1, epsSolver is a technical parameter defining the accuracy 
of a procedure for inverting numerically certain functions - should not be 
modified! The same parameters Xstart and Flavor return the generated 
final value of the x-variable and the final flavor type, while weight is the 
weight of the generated MC event. A detailed history of the evolution is 
recorded inside the object dglap_lo. In particular, the values of all gener- 
ated x- and flavor-variables, including their initial values, are stored there 
(in the m_x:MarkovMC and m_f :MarkovMC matrices). All these parameters 
can be accessed easily with the help of the functions GetParticle and 
GetEventMultiplicity. 



All other evolution types follow exactly the same pattern, as can be seen in 
the demo . cxx file. 



The program demo . cxx prints in the output the final x and the final flavor of 
the first 100 of the generated MC events. After completing the MC generation 
it calculates and prints the average of the MC weight. This average is equal 
to the sum over the final flavors integrated over the final x-variable. This, 
in turn, is almost equivalent to the unitary normalization of the momentum 
distribution functions according to momentum sum rule. It is almost equivalent 
because in the MC program the lower limit e m i n is imposed on the value of the 
generated final a;- variables, x > e min , so there will be a tiny missing piece in 
the sum rule comming from the integral from to e min . Note that this integral 
contains only integrable singularities of the In z-type, so by lowering e m ; n the 
sum rule can be tested to an arbitrary precision. 



On the contrary, the MC solution for x < e m \ n depends on this cut-off and there- 
fore should not be trusted. 

5 In general, Xstart and Flavor will be generated according to some initial parton 
momentum distribution (see Demol) - here they are just set in the code. 
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5.2 More advanced demonstration program in folder Demol 

The library of the classes in the folder MarkovMC is a collection of pure MC 
generators written in a clean and minimalistic way. As shown in the previous 
simple demo, it is easy to use the MC generator objects of these classes. In the 
real life, the pure MC generators are only a small part of a bigger code in which 
they are embedded. Let us call this environmental code the MC application 
(MCAp) and characterize briefly functionality and data structures of such a 
MCAp. The functionality of MCAp typically includes: (a) running many times 
for various input data a MC event generator and storing output data in a form 
of one and more dimensional histograms and/or MC events, in the systems 
with one processor; (b) the same in systems with many processors, many nodes 
(PC farms); (c) visualization and quick analysis of the stored results after the 
MC "production run" is finished, or even while the MC production is still run- 
ning, in particular, (d) comparisons with the stored "benchmark" results from 
the previous "consolidated" versions of the program, (e) comparisons with 
the results of analytical and other non-Monte- Carlo numerical calculations, 
(f) comparisons between the stored MC results from the runs with different 
input data and more. The data structure of MCAp typically features: (A) A 
collection (database) of the input parameters with clear distinction between 
the parameters which are "hardwired" and changed only in very special tech- 
nical tests, the default paremeters which are rarely changed, such that the 
user may normally ignore them, and, finally, the important steering parame- 
ters which are changed often or it is even obligatory to (re-) define them. (B) 
A data base of the output results from the consolidated well-tested versions 
of the program, important results from long CPU time MC runs, outputs 
used in the published works. (C) Some degree of persistency mechanism for 
writing/reading structured data in/from a disk file is absolutely necessary in 
organizing MCAp. This persistency may be limited to data objects, such as 
MC events and histograms, or include a possibility to write into a disk file the 
objects of the random number generator, parts of the MC generator embed- 
ding important member data, or even the complete MC event generator in the 
"ready-to-go" state. 

In the presented distribution package we do not include the full scale MCAp 
environment for the MarkovMC library. However, the demonstration program 
in the folder Demol represents an essential step towards such an infrastructure. 
The demonstration program Demol/MainPr . cxx features to a large extent the 
points (a), (c), (d) and (C) of the above specification list and is meant as an 
useful template for the further development. Let us first overview it briefly, 
and more details will follow later on. 

In this example the Markovian MC generator object m_MMC, being the instance 
of one of the eight classes of the MarkovMC library, is embedded as a member of 
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the object m_MCgen of the container class MMCevol. The object m_MMC is created 
there, filled in with the input data, and used to generate MC events using the 
methods of the MMCevol and MarkovianGen classes. Every MC event generated 
by m_MCgen->m_MMC is made available for histogramming. The object m_MCgen 
is embedded (not created) in the object of the class TRobolA. Methods of this 
class perform several functions: book histograms, transfer input data and a 
pointer of an external random number generator m_RNG object into m_MCgen, 
generate a single MC event, fill-in histograms and write histograms into a disk 
file. The object of the TRobolA class does not contain, however, the main loop 
over the MC event generation, nor creates the mJING and m_MCgen objects. 
The main loop over the MC events is located and managed in a rather special 
way in the main program MainPr, while the objects of the external random 
number generator m_RNG and of the MC generator m_MCgen are created in 
a separate small script Demol/work/Start . C. This scrips, run by ROOT in 
the interpreter mode, defines also all input data of the MC run. Execution of 
Start . C, building, running and stopping the main program (as well as plotting 
the results) is managed by Demol/Makef ile. 

The elements of the above five-level functionality and data structure (Makefile, 
the main program, the TRobolA class, the MMCevol class, the MarkovianGen 
class/library), sketched above, will be described in a more detail in the follow- 
ing. At the first glance, this may look overcomplicated. However, one has to 
remember that the fully functional MCAp, defined in points (a)-(f), (A)-(C) 
above, will unavoidably be rather complicated, especially for running on the 
multinode PC farm. The presented structure exploits years of experience of 
the authors in developing many similar MC infrastructures for running and 
testing MC event generators [34-39], and we hope that it may be useful for 
others as a template to develop it further and/or customize to their own needs. 

Let us now add more details on the Demol source code and its execution. 
All the C++ source code files specific to this demo program are located in 
the folder Demol while the input /output files can be found in its sub folder 
Demol/work. The main() program is located in the file MainPr. cxx. It con- 
sists of three main parts corresponding to the initialization, the generation of 
the MC event series and the final part of the MC program. The same three 
parts are present as three stages in the execution of the program. In the first 
stage, the work/Start. C script creates three objects: (1) an object of the 
MMCevol class, that is the MC generator, (2) an object of an external random 
number generator and (3) an auxiliary object of the class TSemaf . These ob- 
jects are initialized with the default data residing in the class constructors, 
and then they are modified with the run-specific input data contained explic- 
itly in the code of the Start. C script (in particular, a random number seed 
may optionally be redefined at this point). The object of TSemaf contains 
data for administering the main MC loop in the main program. The above 
three objects are written by Start .C into the disk files semaf.root (the object 
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TSemaf ) and mcgen.root (the objects MMCevol and random number generator). 
In the above and in the following the persistency feature of the ROOT system 
is exploited to facilitate writing and reading the entire objects into/from the 
disk files. 

The run-specific input data are all provided in work/Start . C, including all 
steering parameters of the MC generator the user wishes to reset from their 
default values. Examples of such settings together with comments, specific to 
the Demol run, can be found inside work/Start . C. 

After executing Start . C, the main program MainPr comes into action. It 
reads all three stored objects from the disk files semaf.root and mcgen.root 
and creates a new object of the TRobolA class containing all histograms and 
the pointers to the MC generator m_MCgen and the random number genera- 
tor m_RNG. The event generation loop is started. The number of MC events is 
set by the variable NevTot which is read from the TSemaf class object. After 
generating every NGroup events, the partial results of the MC generation are 
stored in a disk (the NGroup variable is also taken from the TSemaf object). 
Before going to the next group of the event generation, the status of the spe- 
cial semaphore-type flag in the TSemaf object in the disk file semaf.root is 
examined. If the status flag is equal to "CONTINUE", then the event generation 
is continued, otherwise, if it is equal to "STOP", then main loop and the MC 
generation is terminated. The latter status flag can be reset by the user in- 
teractively by calling make stop from the subfolder work. At the start of the 
program the status of this flag is always set to "START" and it is redefined to 
"CONTINUE" immediately after generating the first group of the MC events. 
The class called TSemaf is defined in the files (TSemaf . h, TSemaf .cxx). The 
semaphore object contains also some auxiliary information on the generated 
event sample, such as the number of events. Once the event generation loop 
is terminated, the programs enters into the short finalization stage during 
which all the necessary statistics on the MC event sample are calculated and 
all histograms of the TRobolA class are stored in the output file histo.root. 
Histograms are also recorded into histo.root after generating each group of 
MC events. 

The source code of class TRobolA is located in the files (TRobolA . h , TRobolA . cxx) . 
Its main task is to create and fill-in appropriate histograms. This is done in 
three stages: booking, filling and storing of the histograms with the help of 
the class-member functions: HbookO , Production (doublefe) , FileDumpO, 
respectively. More details can be found in the source code. 

The demonstration program Demol /Demo lPr generates a sample of 10 5 weighted 
events for the evolution of the parton momentum distributions from Q m \ n = 
1 GeV to Qmax = 100 GeV. In the course of the MC event generation, the 
program fills-in histograms of the gluon and quark-singlet parton momentum 
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distributions and the corresponding MC weights with the help of ROOT. These 
histograms are stored in the disk file histo.root written in the ROOT format. 
After completing the generation of the series of the MC events, the stored his- 
tograms are plotted using the ROOT graphics programs and compared graphi- 
cally with the histograms pre-computed and stored in the distribution folder. 
The appropriate commands for running this demo program and plotting the 
results are make start and make plot, see the beginning of this section. 

Another C++ program work/draw . cpp is included for analyzing/viewing the 
produced histograms and comparing them with some pre-computed results. 
If the instalation of the MC program is done correctly, the user-produced 
results should agree with the pre-computed ones within statistical errors. This 
program can be run in the interpreter mode of ROOT with the help of a single 
command make plot. 

Let us finally add a few details on embedding the object m_MMC of the MarkovianGen 
class inside the "wrapper" object of the class MMCevol. The object m_MMC 
is created by MMCevol: : Initialize () using the constructor of one of the 
eight classes implementing a specific evolution type. Which one to choose 
is decided by the steering data member m_EvolType of the MMCevol class. 
Then, each MC event is generated by MMCevol: : Generate (), which calls 
m_MMC->GenerateEvent () . Here, the method MarkovianGen: : GenerateEvent () 
does not have any parameters (contrary to the one used in DemoO). This is why 
the getter MarkovianGen : : SetTRange have to be used in MMCevol : : Initialize 
to define the range of the evolution time. Also, prior to m_MMC->GenerateEvent () , 
x and the flavor of the initial parton are generated using the TFoam util- 
ity of ROOT and they are fed into the object m_MMC using the dedicated set- 
ter MarkovianGen: : SetlnitParton. The corresponding part of the code in 
MMCevol: :Generate() looks as follows: 

// Generate primordial/initial parton 
if( m.Foaml != NULL ){ 

m_FoamI->MakeEvent () ; // generate x and parton type 

m_wt = m_FoamI->GetMCwt () ; // get weight 
} else 

m_wt =1.0; 
/// Simulate QCD multiparton evolution 

m_MMC->SetInitParton(m_f lavlni ,m_xlni) ; // Set initial parton 
m_MMC->GenerateEvent() ; // Generate MC event 

m_MMC->GetFinParton(m_Flav,m_X) ; // get final x and parton type 
m_wt *= m_MMC->GetWt() ; // combine MC weight 

The parton momentum distribution of the initial parton is provided to the ob- 
ject m_FoamI of the TFoam class by the dedicated function MMCevol: : Density 
located in Demol/MMCevol . cxx file. The user can easily modify here the shape 
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of the initial distributions. The above organization is quite flexible and allows 
one to introduce into the game more MC generator object^} more MC event 
objects, etc. 



6 Technical tests 



In this section we describe in some detail a variety of tests that we have 
performed on the EvolFMC v. 2 code in order to verify its correctness and 
determine its overall technical precision. To be specific, we have performed 
three different sets of the technical comparisons of the code EvolFMC v. 2: 

(1) with the semianalytical code QCDNuml6 [8], 

(2) with the semianalytical code APCheb40 [12, 13], 

(3) with the previous version of the EvolFMC code: EvolFMC v. 1, and 

(4) between different algorithms within the EvolFMC v. 2. 

We will describe them in turn. The target relative precision of the tests is 
5 x 10~ 4 (half of a per mille). 

At the end of this section we also present the weight distributions for all the 
algorithms and we compare speed of the algorithms. 

As the initial distributions at 1 GeV, for all of the tests we take 

D° G (x) = 1.908 -x- x - 2 (l- xf-\ 
D° sea (x) = 0.6733 ■ x-^tt - x) 7 ' , 



D° u (x) = 2.187 -x-°- 5 (l-x) 3 - , 



D° d ,(x) = 1.230 -x~°- 5 (l -x 



4.0 



val 



The QCD constant A = 0.2457 and Nf = 3. For each of the tests we use 
statistics of the order of 10 10 MC points. 



6.1 Comparison with semianalytical code QCDNuml6 

The previous version of the code, EvolFMC v.l, has been tested against the 
QCDNuml6 code [8] for the standard DGLAP evolution. The demonstrated 
agreement was 5 x 10 -4 for the LO case and 1 x 10~ 3 for the NLO case, 
see [1]. For consistency we have repeated these comparisons for EvolFMC v. 2. 
In Fig. [2] we show the DGLAP NLO evolution up to 10 GeV and 100 GeV 
for both the gluon and quark singlet xDf(x) distributions as well as the ratio 
of the QCDNuml6 to EvolFMC results. The QCDNuml6 results are based on the 

6 Another possibility would be to inherit MMCevol from the MarkovianGen class, 
but such a solution would probably be less flexible. 
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Fig. 2. Left frames : the DGLAP evolutions in the NLO approximation from QCDNuml6 
and EvolFMC (the curves are indistinguishable). Upper curves (magenta and blue): 
the gluon xDg{x) distr. ; lower curves (black and red): the quark xDq(x) distr. 
Right frames: the ratio of QCDNuml6 to EvolFMC for the gluon (magenta) and quark 
(black) distributions. Top frames: the evolution up to 10 GeV. Bottom frames : the evo- 
lution up to 100 GeV. 

extended grid size of 2000 x 600. The agreement is of the order of 5 4- 8 x 10 -4 . 
We will demonstrate in the next subsection that these residual discrepancies 
are to be attributed to QCDNuml6. One has to remember that the comparisons 
with QCDNuml6 are limited to the standard DGLAP evolution only. 



6.2 Comparison with semianalytical code APCheb40 



It is the most independent test of the EvolFMC code. The semianalytical 
method of solving the evolution equations used by APCheb40 is entirely dif- 
ferent from the MC method and also the library of the evolution kernels is 
completely independent. At first, in Fig. |3j we show the comparisons for the 
standard DGLAP NLO evolution (for this comparison we used the older ver- 
sion, 33, of the APCheb code). We show the evolution up to 10 GeV and 100 
GeV for both the gluon and quark singlet xDf(x) distributions as well as the 
ratio of the APCheb33 to EvolFMC results. As we see the agreement is remark- 
able, at the level of 1 x 10~ 4 , limitted by the statistics. This result clarifies the 
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Fig. 3. Left frames : the DGLAP evolutions in the NLO approximation from APCheb33 
and EvolFMC (the curves are indistinguishable). Upper curves (magenta and blue): 
the gluon xDg{x) distr. ; lower curves (black and red): the quark xDq(x) distr. 
Right frames: the ratio of APCheb33 to EvolFMC for the gluon (magenta) and quark 
(black) distributions. Top frames: the evolution up to 10 GeV. Bottom frames : the evo- 
lution up to 100 GeV. 

source of residual disagreement seen in the Fig. [2j 

The advantage of APCheb40 over other semianalytical codes, such as QCDNum [8], 
is that it has the option of the modified-DGLAP evolution built-in, although 
only at the LO level, see [13] for details. For the sake of comparisons we have 
modified the LO kernels in APCheb40 in such a way that the coupling con- 
stant is implemented in the NLO approximation, whereas the ^-dependent 
parts remain in the LO approximation, i.e. the modified kernels are 

X %fp^\ t ^w) = ajVLo(ln( 9 1 "" )+t) 2,p / y ( , ) g 1 _ 2>Ae - t , 

/. , N N (18) 

±x«F™*(t,x,w) = a " Lo(ln( 2 "~ a:)+t) 2,p / y(,)^ >Ae -, 

In Fig. |3]we show the C'-type evolution up to 10 GeV and 100 GeV for both 
the gluon and quark singlet xDf(x) distributions as well as the ratio of the 
APCheb40 to EvolFMC results. The APCheb40 results are based on the interpo- 
lation with 100 Chebyshev polynomials. The agreement is again excellent, in 
most of the x-range below 1 x 10~ 4 . 
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Fig. 4. Left frames : the C'-type evolutions in the modified LO approximation from 
APCheb40 and EvolFMC (the curves are indistinguishable). Upper curves (magenta and 
blue): the gluon xDg(x) distr. ; lower curves (black and red): the quark xDq(x) distr. 
Right frames: the ratio of APCheb40 to EvolFMC for the gluon (magenta) and quark 
(black) distributions. Top frames: the evolution up to 10 GeV. Bottom frames : the 
evolution up to 100 GeV. 

To summarise, the results of comparison with APCheb40 indicate that technical 
precision of EvolFMC is much better than our conservative target of 5 x 10 -4 . 



6.3 Comparison with EvolFMC v.l 



The old EvolFMC v.l has been extensively tested both for the DGLAP (up to 
NLO) and modified-DGLAP (LO only) cases with the overall relative precision 
tag of about 10~ 3 [1, 13]. As an example of the backward compatibility tests 
we show in Fig. [5] the comparison of EvolFMC v.l with EvolFMC v.2 for the 
DGLAP-type NLO evolution and in Fig. [6] for the modified-DGLAP C'-type 
LO evolution, both up to 10 and 100 GeV. As we can see the results agree 
within the statistical errors at the level of 5 x 10 -4 , as desired. 

As we have explained in Introduction, the version v . 1 of EvolFMC has different 
overall structure and organisation of the code as well as different implemen- 
tation of most of the methods. Therefore this comparison can be regarded as 
a comparison with an "almost" independent code (for the limited number of 
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Fig. 5. Left frames : the evolutions from EvolFMC v.l and EvolFMC v. 2 (the curves 
are indistinguishable) for the DGLAP-type NLO evolutions. Upper curves (magenta 
and blue): the gluon xDg{x) distr. ; lower curves (black and red): the quark xDq(x) 
distr. Right frames: the ratio of EvolFMC v. 2 to EvolFMC v.l for the gluon (ma- 
genta) and quark (black) distributions. Top frames: the evolution up to Q = 10 GeV. 
Bottom frames : the evolution up to Q = 100 GeV. 

evolution types, of course). 



6.4 Comparison of different algorithms 



In this section we show the tests of the most advanced evolution: of the C'-type 
in the NLO approximation. We have not found any other code that would solve 
this type of evolution. For this reason we have implemented in the EvolFMC 
code the second, auxiliary, MC algorithm that solves this particular type of 
evolution. The key difference of the auxilliary algorithm with respect to the 
main algorithm is that the entire NLO correction (including modifications in 
the argument of the coupling constant) is introduced as a weight and the 
algorithm itself is based on the LO algorithm described in [28], see [30] for 
details. 

As an example in Fig. [7] we show the comparisons of results of these two NLO 
algorithms of the C'-type for two evolution time limits: 10 GeV and 100 GeV. 
As one can see the agreement is well below the level of 5 x 10 -4 . This result 
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Fig. 6. Left frames : the evolutions from EvolFMC v.l and EvolFMC v.2 (the curves 
are indistinguishable) for modified-DGLAP C'-type LO evolutions. Upper curves (ma- 
genta and blue): the gluon xDq{x) distr. ; lower curves (black and red): the quark 
xDq(x) distr. Right frames: the ratio of EvolFMC v. 2 to EvolFMC v.l for the gluon 
(magenta) and quark (black) distributions. Top frames: the evolution up to Q = 10 
GeV. Bottom frames : the evolution up to Q = 100 GeV. 

we consider as the principal test of the NLO C'-type evolution in the code 
EvolFMC v . 2. Of course, the auxiliary algorithm shares with the previous ones 
a lot of common parts of the code. These common parts however have been 
extensively tested by all the comparisons described in previous subsections. 

We conclude this series of tests with the conservative statement that the 
EvolFMC v.2 code has the overall technical precision of 5 x 10~ 4 ! 



6. 5 Weight distributions and speed of algorithms 



In Fig. [8] we show the weight distributions for all three main algorithms in- 
cluded in the program EvolFMC v.2. The evolutions are of the NLO type 
up to 100 GeV. The weights are well behaved and fall down exponentially. 
The distributions of modified evolutions have shorter tails as compared to 
the DGLAP case. This is the consequence of the finite IR cut-off used in the 
modified kernels, as compared to the infinitesimal one used in the DGLAP 
case. The plots indicate that the conversion to the unweighted events should 
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Fig. 7. Left frames : the modified-DGLAP C'-type evolutions in NLO approximation from 
the two algorithms of EvolFMC v. 2 (the curves are indistinguishable). Upper curves 
(magenta and blue): the gluon xDq{x) distr. ; lower curves (black and red): the quark 
xDq(x) distr. Right frames: the ratio of the two algorithms of EvolFMC v. 2 for the 
gluon (magenta) and quark (black) distributions. Top frames: the evolution up to 10 
GeV. Bottom frames : the evolution up to 100 GeV. 



be quite efficient, especially in the case of moderate target precision level. 
One has to remember that in the current version of the program we did not 
perform any additional weight optimization, as we restricted ourselves to the 
weighted events only. Such an optimization could reduce the tails of the weight 
distributions even further. 



Finally, in Fig. [9] we compare speed of all five algorithms included in the 
program EvolFMC v. 2 both in LO and NLO cases. As we see, despite their 
complexity and the presence of internal one-dimensional numerical integration, 
the NLO algorithms are slower only by a factor of few as compared to the very 
fast LO ones. The slowest one is the B'-type algorithm. Let us note that the 
algorithms are only partly optimized with respect to speed. 
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Fig. 8. Weight distributions for NLO evolution to 100 GeV (normalized to one). 
Left frames : gluon. Right frames: quark. Top frames: DGLAP evolution. Central frames : 
B'-type evolution. Bottom frames : C'-type evolution. 

7 Summary 



In this paper we have presented the program EvolFMC v . 2. It is a Monte Carlo 
generator that solves some of the QCD evolution equations for the parton mo- 
mentum distributions in the weighted event mode. In the current version we 
have implemented the DGLAP evolution in the LO and NLO approximations 
and two modified-DGLAP type evolutions. These modifications include the 
replacement of the argument of the coupling constant by a more complicated 
functions of the x and t variables. In one case (C) this function reconstructs 
the transverse momentum of the emitted partons and in the other case (B') 
the approximate transverse momentum. Both of these modified-DGLAP evo- 
lutions are implemented in the LO as well as the NLO approximation. The 
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Fig. 9. Comparison of the speed of all five algorithms. Plotted times (in seconds) refer 
to generation of 10 6 events and evolution to 100 GeV. First three blocks illustrate the 
LO evolution (DGLAP, B' and C, respectively), following five blocks the NLO case 
(DGLAP, B'-auxilary, B'-main, C'-auxiliary and C'-main, respectively). 

code is written in the C++ language and has the modular structure, easy to 
extend to additional evolution types. The main limitation of the code is the 
fact that quarks are treated as massless. We have studied the technical pre- 
cision of the program by means of extensive comparisons with the non-MC 
numerical program APCheb40, with the previous version of the code, EvolFMC 
v.l, and by comparing two independent algorithms for the same evolution 
implemented in EvolFMC v . 2. These tests have proved that the technical pre- 
cision of the EvolFMC v. 2 is at least 5 x 10~ 4 . To our knowledge it is the first 
so-precise solution of the QCD evolution equations by means of the Monte 
Carlo methods. 
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9 Appendix 

In this Appendix we list the source code of the simple demo . cxx program from 
the DemoO folder. 
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#include<stdio . h> 
#include<math . h> 



#include "rndm.h" 

#include "bprim_nlo.h" 
#include "bprim_lo .h" 
#include "bprim_nlo_aux .h" 
#include "cprim_lo .h" 
#include "cprim_nlo.h" 
#include "cprim_nlo_aux.h" 
#include "dglap_lo.h" 
#include "dglap_nlo.h" 



#def ine 


DGLAP. 


LO 







#def ine 


DGLAP. 


NLO 




1 


#def ine 


CPRIM. 


.LO 




2 


#def ine 


CPRIM. 


.NLO 




3 


#def ine 


CPRIM. 


NLO. 


_AUX 


4 


#def ine 


BPRIM. 


.LO 




5 


#def ine 


BPRIM. 


.NLO 




6 


#def ine 


BPRIM. 


.NLO. 


.AUX 
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double Lambda = 0.25; 

double numberOf Flavor = 3; 
double Qmin = 1 ; 

double Qmax = 100; 

int TypeOf Generator = 0; 

rndm * RNgen; 

int numberOf Events = 100; 

int main(int argc, char *argv[]) 
{ 

if (argc>l) 

TypeOf Generator = (atoi (* (++argv) ) <8) ? atoiOargv) : 0; 

RNgen = new rndm() ; 
double TO = log(Qmin) ; 
dglap_lo * dglap_ll 

= new dglap_lo( TO, Lambda, number Of Flavor , RNgen) ; 
dglap_nlo * dglap_nll 

= new dglap_nlo( TO, Lambda, number Of Flavor, RNgen) ; 
CPrim_lo * cprim_ll 

= new CPrim_lo( TO, Lambda, number Of Flavor, RNgen) ; 
CPrim_nlo * cprim_nll 
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= new CPrim_nlo( TO, Lambda, number Of Flavor, RNgen) ; 
CPrim_nlo_aux * cprim_nll_aux 

= new CPrim_nlo_aux (TO, Lambda, numberOf Flavor ,RNgen) ; 
BPrim_nlo * bprim_nlo 

= new BPrim_nlo( TO, Lambda, number Of Flavor, RNgen) ; 
BPrim_lo * bprim_lo 

= new BPrim_lo( TO, Lambda, number Of Flavor , RNgen) ; 
BPrim_nlo_aux * bprim_nlo_aux 

= new BPrim_nlo_aux (TO, Lambda, numberOf Flavor , RNgen) ; 

dglap_ll->printLogo() ; 

int init_Flavor = 0; 
double init_Xstart = 1.0; 

printf ( " *******************Starting parameters* ****** **********\n" ) ; 
printf ("Qmin=°/4.2f\n", Qmin) ; 
printf ("Qmax= > / 4.2f\n" , Qmax) ; 

printf ("Type of Equation: °/ d\n" , TypeOf Generator) ; 
printf ("QCDLambda=%2.4f\n" , Lambda) ; 
printf ( "numberOf Flavor s=%d\n" , numberOf Flavor) ; 
printf ("Events to generate : °/ d\n" , numberOf Event s) ; 

printf ( " *******************Initial Conditions* ******* **********\n" ) ; 
printf ("Flavor : %d\n" , init_Flavor) ; 
printf ("Xstart: °/ 1.2f\n", init_Xstart) ; 

printf ( " *************************Events************************\n" ) ; 

double TStop = log(Qmax) ; 
double Tc; 

double sumwt = 0.0; 
double sumwt2= 0.0; 
int NoEvents =0; 

for (int i=0; KnumberOf Events ; i++) 
{ 

int Flavor = init_Flavor; 
double Xstart = init_Xstart ; 

double weight = 1.0; 

switch (TypeOf Generator) 
{ 

case CPRIM_NL0: 

cprim_nll->GenerateEvent (Tc , Xstart, weight, TStop, Flavor); 
break; 
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case CPRIM_NLD_AUX: 

cprim_nll_aux->GenerateEvent (Tc , Xstart, weight, TStop, Flavor); 

break; 
case CPRIM_L0: 

cprim_ll->GenerateEvent (Tc , Xstart, weight, TStop, Flavor); 
break; 
case DGLAP_L0: 

dglap_ll->GenerateEvent (Tc , Xstart, weight, TStop, Flavor); 
break; 
case DGLAP_NL0: 

dglap_nll->GenerateEvent (Tc , Xstart, weight, TStop, Flavor); 
break; 
case BPRIM_L0: 

bprim_lo->GenerateEvent (Tc , Xstart, weight, TStop, Flavor); 
break; 
case BPRIM_NL0 : 

bprim_nlo->GenerateEvent (Tc , Xstart, weight, TStop, Flavor); 
break; 
case BPRIM_NLD_AUX: 

bprim_nlo_aux->GenerateEvent (Tc , Xstart, weight, TStop, Flavor); 
break; 

} 

sumwt =sumwt + weight; 
sumwt2=sumwt2 + weight*weight ; 
NoEvents = NoEvents+1; 
if (NoEvents<100) 
{ 

printf ("X=°/„1.6f , weight=%l . 6f , f lavor=7„d\n" , Xstart, weight, Flavor); 

} 

} 

double integral=sumwt/NoEvents ; 
double error = sqrt( (sumwt2/NoEvents 

- (sumwt /NoEvents) * (sumwt /NoEvents) ) /NoEvents); 

printf ( " ***********************0utput**************************\n" ) ; 
printf ("\n [sum_f =f lavor] [int_eps"l dx] D_f(x,t)= %1.6f +- %1.6f\n\n" 

.integral, error); 

return 0; 
} 
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